{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###### Content under Creative Commons Attribution license CC-BY 4.0, code under BSD 3-Clause License © 2018 parts of this notebook are from [Derivative Approximation by Finite Differences](https://www.geometrictools.com/Documentation/FiniteDifferences.pdf) by David Eberly,  additional text and SymPy examples by D. Koehn, notebook style sheet by L.A. Barba, N.C. Clementi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Execute this cell to load the notebook's style sheet, then ignore it\n",
    "from IPython.core.display import HTML\n",
    "css_file = '../../style/custom.css'\n",
    "HTML(open(css_file, \"r\").read())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Generalization of Taylor FD operators\n",
    "\n",
    "In the last lesson, we learned how to derive a high order FD approximation for the second derivative using Taylor series expansion. In the next step we derive a general equation to compute FD operators, where I use a detailed derivation based on [\"Derivative Approximation by Finite Differences\" by David Eberly](https://www.geometrictools.com/Documentation/FiniteDifferences.pdf)\n",
    "\n",
    "## Estimation of arbitrary FD operators by Taylor series expansion\n",
    "\n",
    "We can approximate the $d-th$ order derivative of a function $f(x)$ with an order of error $p>0$ by a general finite-difference approximation:\n",
    "\n",
    "\\begin{equation}\n",
    "\\frac{h^d}{d!}f^{(d)}(x) = \\sum_{i=i_{min}}^{i_{max}} C_i f(x+ih) + \\cal{O}(h^{d+p})\n",
    "\\end{equation}\n",
    "\n",
    "where h is an equidistant grid point distance. By choosing the extreme indices $i_{min}$ and $i_{max}$, you can define forward, backward or central operators. The accuracy of the FD operator is defined by it's length and therefore also the number of \n",
    "weighting coefficients $C_i$ incorporated in the approximation. $\\cal{O}(h^{d+p})$ terms are negelected. \n",
    "\n",
    "Formally, we can approximate $f(x+ih)$ by a Taylor series expansion:\n",
    "\n",
    "\\begin{equation}\n",
    "f(x+ih) = \\sum_{n=0}^{\\infty} i^n \\frac{h^n}{n!}f^{(n)}(x)\\nonumber\n",
    "\\end{equation}\n",
    "\n",
    "Inserting into eq.(1) yields\n",
    "\n",
    "\\begin{align}\n",
    "\\frac{h^d}{d!}f^{(d)}(x) &= \\sum_{i=i_{min}}^{i_{max}} C_i \\sum_{n=0}^{\\infty} i^n \\frac{h^n}{n!}f^{(n)}(x) + \\cal{O}(h^{d+p})\\nonumber\\\\\n",
    "\\end{align}\n",
    "\n",
    "We can move the second sum on the RHS to the front\n",
    "\n",
    "\\begin{align}\n",
    "\\frac{h^d}{d!}f^{(d)}(x) &= \\sum_{n=0}^{\\infty} \\left(\\sum_{i=i_{min}}^{i_{max}} i^n C_i\\right) \\frac{h^n}{n!}f^{(n)}(x) + \\cal{O}(h^{d+p})\\nonumber\\\\\n",
    "\\end{align}\n",
    "\n",
    "In the FD approximation we only expand the Taylor series up to the term $n=(d+p)-1$\n",
    "\n",
    "\\begin{align}\n",
    "\\frac{h^d}{d!}f^{(d)}(x) &= \\sum_{n=0}^{(d+p)-1} \\left(\\sum_{i=i_{min}}^{i_{max}} i^n C_i\\right) \\frac{h^n}{n!}f^{(n)}(x) + \\cal{O}(h^{d+p})\\nonumber\\\\\n",
    "\\end{align}\n",
    "\n",
    "and neglect the $\\cal{O}(h^{d+p})$ terms\n",
    "\n",
    "\\begin{align}\n",
    "\\frac{h^d}{d!}f^{(d)}(x) &= \\sum_{n=0}^{(d+p)-1} \\left(\\sum_{i=i_{min}}^{i_{max}} i^n C_i\\right) \\frac{h^n}{n!}f^{(n)}(x)\\\\\n",
    "\\end{align}\n",
    "\n",
    "Multiplying by $\\frac{d!}{h^d}$ leads to the desired approximation for the $d-th$ derivative of the function f(x):\n",
    "\n",
    "\\begin{align}\n",
    "f^{(d)}(x) &= \\frac{d!}{h^d}\\sum_{n=0}^{(d+p)-1} \\left(\\sum_{i=i_{min}}^{i_{max}} i^n C_i\\right) \\frac{h^n}{n!}f^{(n)}(x)\\\\\n",
    "\\end{align}\n",
    "\n",
    "Treating the approximation in eq.(2) as an equality, the only term in the sum on the right-hand side of the approximation that contains $\\frac{h^d}{d!}f^{d}(x)$ occurs when $n = d$, so the coefficient of that term must be 1. The other terms must vanish for there to be equality, so the coefficients of those terms must be 0; therefore, it is necessary that\n",
    "\n",
    "\\begin{equation}\n",
    "\\sum_{i=i_{min}}^{i_{max}} i^n C_i=\n",
    "\\begin{cases}\n",
    "0, ~~ 0 \\le n \\le (d+p)-1 ~ \\text{and} ~ n \\ne d\\\\\n",
    "1, ~~ n = d\n",
    "\\end{cases}\\nonumber\\\\\n",
    "\\end{equation}\n",
    "\n",
    "This is a set of $d + p$ linear equations in $i_{max} − i_{min} + 1$ unknowns. If we constrain the number of unknowns to be $d+p$, the linear system has a unique solution. \n",
    "\n",
    "- A **forward difference approximation** occurs if we set $i_{min} = 0$\n",
    "and $i_{max} = d + p − 1$. \n",
    "\n",
    "- A **backward difference approximation** can be implemented by setting $i_{max} = 0$ and $i_{min} = −(d + p − 1)$.\n",
    "\n",
    "- A **centered difference approximation** occurs if we set $i_{max} = −i_{min} = (d + p − 1)/2$ where it appears that $d + p$ is necessarily an odd number. As it turns out, $p$ can be chosen to be even regardless of the parity of $d$ and $i_{max} = (d + p − 1)/2$.\n",
    "\n",
    "We could either implement the resulting linear system as matrix equation as in the previous lesson, or simply use a `SymPy` function which gives us the FD operators right away."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import SymPy libraries\n",
    "from sympy import symbols, differentiate_finite, Function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define symbols\n",
    "x, h = symbols('x h')\n",
    "f = Function('f')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Actually, the underlying algorithm also supports variable grid spacings, because it is not based on Taylor series expansion but Lagrange polynomials. For more details, I refer to the paper [\"Calculation of weights in finite difference formulas\" by Bengt Fornberg](https://amath.colorado.edu/faculty/fornberg/Docs/sirev_cl.pdf)\n",
    "\n",
    "An alternative to using `SymPy` for the calculation of FD operator weights is the **DEVITO** package:\n",
    "\n",
    "[\"https://www.devitoproject.org/\"](https://www.devitoproject.org/)\n",
    "\n",
    "It does not only calculate FD stencils, but also automatically optimizes the performance for the given hardware."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## What we learned:\n",
    "\n",
    "* How to compute Finite-Difference operators of arbritary derivative and error order\n",
    "* Symbolic computation of FD operators with `SymPy`"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
